/*  README:
*       Algoritmo de Gauss Legendre(Paralelo)
*
*           Algoritmo que converge na ordem de 2 casas decimais com precisao de 10.000.000 casas decimais,
*           portanto temos 25 iteracoes para seu calculo na main;
*           para alterar a quantidade de iteracoes deve-se alterar o iterador de 25 para o desejado na funcao;
*           o tempo eh fornecido em milisegundos,; defina a precisao na variavel "prec"
*
*/


#include<stdio.h>
#include<stdlib.h>
#include<gmp.h>
#include<pthread.h>
#include<time.h>

typedef struct{
    mpf_t a0, a1;
    mpf_t b0, b1;
    mpf_t t0, t1;
    mpf_t p0, p1;
    mpf_t pi;
}gmpf_t;

//Inicializa a0 com 1
void *inicializaA(mpf_t *a0){
    mpf_set_str(*a0, "1",10);
}

//Inicializa b0 com 1/sqrt(2)
void *inicializaB(mpf_t *b0){
    mpf_t aux;
    mpf_init(aux);

    mpf_set_str(*b0, "1",10);
    mpf_set_str(aux, "2",10);

    mpf_sqrt(aux,aux);

    mpf_div(*b0, *b0, aux);

    mpf_clear(aux);
}

//Inicializa t0 com 1/4
void *inicializaT(mpf_t *t0){
    mpf_t aux;
    mpf_init(aux);

    mpf_set_str(*t0, "1",10);
    mpf_set_str(aux, "4",10);

    mpf_div(*t0, *t0, aux);

    mpf_clear(aux);
}

//Inicializa p0 com 1
void *inicializaP(mpf_t *p0){
    mpf_set_str(*p0, "1",10);
}

void *proximoA(gmpf_t *vars){
    mpf_t aux;
    mpf_t aux2;
    mpf_init(aux);
    mpf_init(aux2);

    mpf_add(aux2, vars->a1, vars->b1);
    mpf_set_str(aux, "2",10);

    mpf_div(vars->a1, aux2, aux);

    mpf_clear(aux);
    mpf_clear(aux2);
}

void *proximoB(gmpf_t *vars){
    mpf_t aux;
    mpf_init(aux);

    mpf_mul(aux, vars->a0, vars->b0);
    mpf_sqrt(vars->b1, aux);

    mpf_clear(aux);
}

void *proximoT(gmpf_t *vars){
    mpf_t aux;
    mpf_init(aux);

    mpf_sub(aux, vars->a1, vars->a0);
    mpf_pow_ui(aux, aux, 2);
    mpf_mul(aux, aux, vars->p0);

    mpf_sub(vars->t1, vars->t0, aux);

    mpf_clear(aux);
}

void *proximoP(gmpf_t *vars){
    mpf_t aux;
    mpf_init(aux);

    mpf_set_str(aux, "2", 10);

    mpf_mul(vars->p1, vars->p0, aux);
}

void calculaPi(gmpf_t *vars){
    mpf_t aux;
    mpf_t aux1;
    mpf_t aux2;
    mpf_init(aux);
    mpf_init(aux1);
    mpf_init(aux2);

    mpf_add(aux, vars->a0, vars->b0);
    mpf_pow_ui(vars->pi, aux, 2);

    mpf_set(aux1, vars->t0);
    mpf_set_str(aux2, "4", 10);
    mpf_mul(aux, aux1, aux2);

    mpf_div(vars->pi, vars->pi, aux);
}

int main(){
    clock_t tInicio, tFim;
    tInicio = clock();

    int i, prec = 33211928;//10.000.000 * log(10)/log(2)
    pthread_t  thread[9];
    gmpf_t var;

    mpf_set_default_prec(prec);

    mpf_init(var.a0);
    mpf_init(var.b0);
    mpf_init(var.t0);
    mpf_init(var.p0);
    mpf_init(var.a1);
    mpf_init(var.b1);
    mpf_init(var.t1);
    mpf_init(var.p1);
    mpf_init(var.pi);

    pthread_create( &thread[0], NULL, (void *)inicializaA, (void *)&(var.a1)); /*Inicializa a0*/
    pthread_create( &thread[1], NULL, (void *)inicializaB, (void *)&(var.b1)); /*Inicializa b0*/
    pthread_create( &thread[2], NULL, (void *)inicializaT, (void *)&(var.t1)); /*Inicializa t0*/
    pthread_create( &thread[3], NULL, (void *)inicializaP, (void *)&(var.p1)); /*Inicializa p0*/

    pthread_join( thread[0], NULL);
    pthread_join( thread[1], NULL);
    pthread_join( thread[2], NULL);
    pthread_join( thread[3], NULL);

    for(i=0; i<25; i++){
        mpf_init_set(var.a0, var.a1);
        mpf_init_set(var.b0, var.b1);
        mpf_init_set(var.t0, var.t1);
        mpf_init_set(var.p0, var.p1);

        pthread_create( &thread[4], NULL, (void *)proximoA, (void *)&(var)); /*Inicializa a0*/
        pthread_create( &thread[5], NULL, (void *)proximoB, (void *)&(var)); /*Inicializa b0*/
        pthread_join( thread[4], NULL);
        pthread_join( thread[5], NULL);

        pthread_create( &thread[6], NULL, (void *)proximoT, (void *)&(var)); /*Inicializa t0*/
        pthread_create( &thread[7], NULL, (void *)proximoP, (void *)&(var)); /*Inicializa p0*/
        pthread_create( &thread[8], NULL, (void *)calculaPi, (void *)&(var)); /*Inicializa p0*/

        pthread_join( thread[6], NULL);
        pthread_join( thread[7], NULL);
        pthread_join( thread[8], NULL);

        gmp_printf("Interacao %d: PI = %.30Ff \n",i+1, var.pi);
    }

    tFim = clock();
    printf("\nNumero Total de Interacoes: %d    Tempo Decorrido:%ld ms\n\n", i,(tFim-tInicio)/(CLOCKS_PER_SEC/1000));


    return 0;
}

